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(— ! , Abstract. We introduce a new approach to assess the error of control problems we 

f*> i' aim to optimize. The method offers a strategy to define new control pulses that are 

not necessarily optimal but still able to yield an error not larger than some fixed a 
priori threshold, and therefore provide control pulses that might be more amenable 
r^ ' for an experimental implementation. The formalism is applied to an exactly solvable 

O^. model and to the Landau-Zener model, whose optimal control problem is solvable only 

numerically. The presented method is of importance for applications where a high 
04 ■ degree of controllability of the dynamics of closed quantum systems is required. 
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1. Introduction 

The theory of optimal control (OC) that has been mathematically formulated in the last 
century by the seminal works of Pontryagin, Bellman, Kalman, Stratonovich [H [21 13] 
has been instrumental for the achievement of highly reliable electronic devices used to 
control, for instance, mechanical systems such as airplanes, cars, etc., but also to control 
chemical reactions or to design ultra-fast laser pulses for manipulating molecules (e.g., 
to break a certain bond while leaving other bonds intact [1]), and today even to optimize 
(stochastic) financial analyses [5]. 

The topic has recently attracted the attention of physicists working in quantum 
information and computation science, because of the need to engineer accurate protocols. 
To this aim different numerical techniques have been devised in order to minimize (or 
maximize) some performance criterion or, alternatively called, objective functional. 
We mention the most used ones in open-loop quantum control: the Krotov iterative 
method P, El [7] and the gradient ascent pulse engineering algorithm |H]. Although 
these are powerful tools for the search of OC pulses they do not provide an assessment 
of the tolerable error against distortions and do not guarantee to obtain control 
pulses easily realizable in the laboratory. Such an issue is of paramount importance 
for quantum information processing (QIP), since the error allowed by fault-tolerant 
quantum computation ranges between 0.01% to fractions of a percent P, fTO] . 

A possible (empirical) approach relies on applying an arbitrary distortion to the 
OC solution and then looking at the error that it produces on the objective functional, 
or by selecting a region in the Hilbert space that is robust against noise (decoherence 
free-subspace) [11], whose existence follows from the symmetry properties of the noise. 
Recently, a more systematic methodology based on the Hessian analysis of the cost 
functional has been proposed [12] or by using an improved genetic algorithm in the 
presence of control noise [13] . 

The aim of this work is to provide an alternative method to evaluate to which 
extent a given OC scheme can tolerate errors. The method is based on the Hessian 
approximation of the cost functional, as in Ref . [12] , but it is applicable to any system 
Hamiltonian H{ut) and Hilbert space dimension. Such arbitrariness is important in 
several circumstances where either the control pulse Ut = u{t) does not appear linearly 
in H{ut) or where the state to be controlled is an auxiliary state (e.g., the motional 
state of an atom [H]) and not the quantum bit itself (e.g., an atomic internal state). 
The control of such states is relevant not only for several QIP implementations, but also 
for quantum metrological purposes, where the control of large quantum superpositions 
may increase the sensitivity of precise measurements. 

Even though usually it is not possible to know analytically the OC pulse, we 
underscore that our assumption is that the parameter obtained with some numerical 
algorithm is very close to the global optimum. More precisely, the error on the cost 
functional obtained with the numerically found OC pulse has to be much smaller than 
the error allowed by the process we are interested to optimize. Besides the interest on 




Figure 1. (Color online) Pictorial representation of the optimal trajectory |?/;°) (red- 
thick line) obtained with the OC pulse u° and trajectories (thin lines) for non-optimized 
control pulses. The shaded (green) area represents the portion of Hilbert space within 
which the cost functional i < J (see also text), that is, the subset of state vectors 
close to the goal state IV'g)- 



its own, we believe that our approach might be of importance for experiments, where, 
typically, optimal pulses are extremely difficult to achieve. To this aim, our method 
could help to find easily implementable control signals (EICS), while still being able 
to satisfactorily fulfill the performance criterion we are interested in. Here with "easily 
implementable control signals" reference is made to pulses that can be utilized to control 
an experiment at the quantum level. More precisely, since nowadays the experiments 
are typically controlled by a computer, an obvious requirement for the control signal is 
that its Fourier spectrum has to match the bandwidth of the transducer or it can not 
vary faster than the clock frequency of the processor. Besides this, since the computer 
during the course of the experiment controls some device (e.g., the applied voltage on 
electrodes or electric current [HI [15]) the control pulse has, for instance, to take into 
account the bandwidth of those devices. These conditions might be not satisfied by the 
optimal control pulse obtained with the aforementioned optimization algorithms. Even 
though, recently, some extensions of those optimization methods in order to include 
spectral constraints on the control pulses have been made p^fTT]. these are not always 
easy to be handled, especially when the dynamics of a many-body quantum system is 
concerned. 

2. The method 

Let us briefiy review what is the purpose and what are the methodologies adopted 
commonly in the quantum control research area concerning closed quantum systems 
and the link between OC theory and analytical mechanics. The usual problem is the 
engineering of a system Hamiltonian H{ut) such that the initial state \ipir^ at time t = to 
of the quantum system under consideration is brought at time t = T to some desired 
goal state \ipg) (see also Fig. [1]). To this aim there are two (equivalent) techniques 
for searching optimal pulses at our disposal: the variational method and dynamical 
programming. In both cases the goal is the minimization of the cost functional 

J[to,ut,iJt] = G(^t) + / dtC{t,ut,iJt) (1) 
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over all admissible control pulses Ut and state trajectories \ipt) = I'ipit)) given a certain 
initial state iV'in)- Here by admissible we mean any Ut for which the Schrodinger equation 
is well-defined and has a unique solution \ijjt) given the initial condition \ipin)- The first 
term on the r.h.s. of Eq. ([1]) is the terminal functional, e.g., the overlap infidelit 
1 — \{ipg\il)T)\'^- The second term provides additional constraints on the control pulse 
In the variational method the additional constraint Re{/^ dt[{xt\'(pt) + ^iXtlH {ut)\ipt)]} 
is introduced [18], where we set h = 1 and \xt) is a Lagrange multiplier often referred 
to as costate, which ensures that the state lipt) satisfies the Schrodinger equation. The 
search of an extremum of the cost functional produces a set of equations for the state, 
costate and the control pulse. 

On the other hand, dynamical programming, based on the Bellman's optimality 
principle [19], produces an equation, the so called Hamilton- Jacobi-Bellman equation, 
for the optimal cost S(t,4't) '■= infu^ J[t, Wt,'?/'*]- This equation is expressed in formally 
the same way as the Hamilton- Jacobi equation of classical mechanics, but with the 
difference that it is propagated backwards in time. The role of the Hamilton function in 
classical mechanics is played in control theory by the (quantum) Pontryagin Hamiltonian 
H(x, "^A) := snp^^{Re[i{xt\H{ut)\^t)] — ^i't,ut,ipt)} [201 EI], which is not to be confused 
with the system Hamiltonian H{ut) that we control through Ut- From the Hamilton- 
Jacobi-Bellman equation it is possible to retrieve the equations of motion for the 
state and the costate, the same as in the variational approach, which look, given 
the aforementioned Pontryagin Hamiltonian, formally as the Hamilton equations for 
the phase space variables {q,p) in analytical mechanics. Beside this, the solution of 
the Hamilton- Jacobi-Bellman equation, defined through the solution of the latter (the 
Hamilton boundary- value problem), is given by |21j : 

S(to,V'o) = G(^t) + / dt[{ijt\xt) - H(xt,^i)], (2) 

J to 

which is very similar to the action in classical mechanics. We note, however, that the 
variational approach, based on the so-called Pontryagin maximum principle p3], yields 
necessary and sufficient conditions for local minima, whereas dynamical programming 
produces results that are globally optimal. 

Motivated by this analogy between OC theory and analytical mechanics we can 
view the cost functional as an "action functional" . Indeed, in analogy to the Hamilton's 
principle where the actual evolution of a classical system is an extremum of the action 
functional, which produces the well-known Euler-Lagrange equations of motion, an 
extremum of J produces the OC equations for the state and the costate. 

Now, let us describe how our method works. To begin with, we fix the desired cost, 
J', that the system (at least) has to attain. Since the advantage of optimizing quantum 
dynamics is connected with the possibility of reaching \i/jg) by exploiting the interference 
of several paths in the space of control parameters U, we define a path integral in such 
a space. To this aim, we introduce the weight K = J 'D[ut]e"t 1*"''''*^ with 1)[ut\ being 

I For instance, the "laser electric field fluence" with C(uj) = Uj , where Ut is an electric field amplitude. 



some suitable measure on the space U. The form of this weight resembles the Feynman 
propagator, which is motivated by the previous discussed analogy between analytical 
mechanics and OC theory, and by the expression ([2]) for the optimal cost. The choice 
of the exponential, however, does not emerge from fundamental physical requirements. 
We think that any well-behaved function peaked around the optimal control pulse will 
allow to estimate the robustness of an optimal control problem. This conjecture is based 
on the observation, see the discussion in the following, that according to our analysis 
the curvature of the cost functional around the optimal solution is the relevant quantity 
to analyze. 

The numerical factor at given in the weight K, which we shall refer to "infidelity 
tolerance" , has the following property: when J approaches (from above) the minimum 
of J, then ttt — > 0. Hence, at is the analogous of h in Feynman path integral. Indeed, 
when at — !■ we retrieve the OC equations for the state and the costate we mentioned 
before. Besides this, we assume that at exists and it is unique for a given OC problem, 
regardless of the form of the distortion 5ut = Ut — u°. The uniqueness is first of all a 
necessary condition or else our method would be ineffective. There is, however, a more 
deep reason why we make such an assumption. This is more easily understood in the 
case where C = in Eq. ([1]). In such a scenario the cost functional relies solely upon the 
state at the final time T, that is, the cost functional will depend on some time integral of 
Ut over [to, T] (and eventually on other time integrals for its time derivatives over [to, T] 
depending on the particular control problem one is interested in). Thus, minimizing the 
cost functional means to find the right conditions for those integral functionals which 
depend only on time T. Hence, the form of Ut in the interval [to,T) does not matter, if 
those functionals of the control pulse and its time derivatives fulfill the right conditions 
at the final time T. In this respect the infidelity tolerance has to be independent from 
the distortion we utilize, that is, it is unique. This is also what is shown by the analysis 
carried out on the examples we will discuss later in the paper. 

Even though it relies on the particular control problem we have at hand, from the 
above outlined discussion, at least in the most relevant cases for QIP where C = in 
Eq. ([1]), the set of EICS, A, is dense, since what matters is the fulfillment of the right 
conditions at the final time T. For instance, in the first example we are going to consider 
in the next section, that is, the optimal transport of a particle confined in a moving 
harmonic trap, Ref. p2] has showed that if u° is the optimal solution then u° = u° + au° 
Va G M is optimal as well (here JT" = 0, see also Sec. [3]). Thus, there is a group of 
optimal solutions, parametrized by the continuous variable a, which is topologically 
dense. Similarly, this will occur for JT 7^ 0, whose set of control pulses forms another 
(dense) subset it in U. Each of these possible control signals will have a precise spectrum 
that has to be within the bandwidth of it, namely the largest bandwidth of the elements 
of ^. Now, if the EICS needs to have a specific bandwidth, then one has to select from 
it the Ut that have the right bandwidth, namely restrict (continuously) the bandwidth 
of il such that the right subset of 11 becomes A^ even though such a procedure might 
produce an empty set. 
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To be specific let us set, for tlie sake of simplicity, to = and consider only one 
control pulse Ut : [0, T] — )■ M, that is, only one control pulse is applied to the system we 
aim to steer during the time interval [0, T]. We also assume that the state \ipt) obeys the 
Schrodinger equation with the time-dependent Hamiltonian operator H{ut) and initial 
condition \iPq) = \x(:{0)). Contrarily to ([T]), where the state ip is an independent variable, 
hereafter we render explicit the dependence of tp on the control parameter ut- Thus, 
we introduce the reduced cost functional J'[n] := J[n,^(n)]. Then by performing the 
Taylor expansion of J'[n] around u° to second order in 6ut we obtain 

J'[u] ~ J'[u°] + ^6uH6u^ + 0{\\ 6u f ), (3) 

where || ■ || is some norm in U, J'[u°] is the minimum, which, without loss of generality, 
we will set to zero, 6vJ is the transposed of the vector 6u = {6u2, ■ ■ ■ , 6u]\f^i), and the 
time interval [0,T] is divided in iV — 1 equal parts At = T/{N — 1). The Hessian H 
is a real, symmetric, and positive defined matrix, and therefore it can be diagonalized. 
Additionally, we assume the boundary conditions Suq = Sut = 0, that is, the values 
uq and Ut are fixed. We underscore that henceforth we shall work with the discretized 
system time evolution since in most cases the optimization of a given control problem 
is performed numerically, and therefore it is discretized, but most importantly, because 
any experiment is performed with a finite number of control time steps. 

Given that, let us define a suitable norm in U for 6ut by means of K, in order 
to obtain a quantitative appraisal of the tolerated error by a given OC scheme. A 
suitable choice for such a norm is given by: (J') = | JS)[Mt]J'[Mt] exp(^J'[ut])| with 
f T)[ut] = \im.]\[^oo f duN-i ■ ■ ■ J du2- Because of the second order approximation in 
Eq. (|3]), the integrals appearing in (J') over the control pulse at different times can be 
replaced (by making a change of integration variables) with Su^ V/c = 2, . . . ,N — 1. 
This norm is indeed an "average" in U of the reduced cost functional itself with the 
exponential function being a "probability distribution". A similar approach has been 
introduced by Rabitz [23] where the terminal functional G in Eq. ([1]) is averaged over 
a distribution function P{6ut). Within our method P{6ut) can be identified with the 
exponential function in (J'). We remark that we are not interested in some specific 
noise model, but rather to identify the portion ii gU which enable us to satisfy J' < JT". 
Basically, il is determined by (J'), which only relies on u° and J^, and any kind of 
noise has to yield pulses such that Ut G il. Hence, our approach is applicable to any 
OC problem and cost functional ((H). Furthermore, we note that in the next we shall 
consider only real- valued Ut (e.g., an electric current [H]), but we underscore that the 
path integral can be easily generalized to complex-valued Ut like the amplitude and the 
phase of a laser field. 

Close to the optimal solution we can approximate J' with its second order expansion, 
and therefore the norm becomes 




(j')-A/7^"?y^^> (4) 



where N'i7'^ = \ Ylj^tk I ^^j^^ j?j/(2"t)| ^^^^ normalization factors, and A^ are the non-zero 
eigenvalues of H with M < N — 2. The above formula makes good sense, because when 
at -^ also (J') -)■ 0. 

In order to apply the above outlined formalism to some concrete example we shall 
consider hereafter y[ut] = 1 — ^H^t), with f{ipT) = K'^gl'^r)!^- We note that the state 
lipx) is implicitly depending on the whole history of Ut Vt G [0,T). Assuming that 
F{ut,ipt{ut)) is a differentiable functional of its arguments we have 



J'K] = 1 - F(^o) + 2 / dtlm[{^,\^t){^t\H{ut)\iJ,)]. (5) 

Jo 

Here we used the fact that f{ipT) = Fii'o) + Jq dt^('?/^f). The most difficult part of 
our method is the computation of the Hessian matrix H. To this aim we have two 
possibilities at our disposal: either we estimate H by means of the Broyden-Fletcher- 
Goldfarb-Shanno (BFGS) formula [24] or we compute the first and second derivatives of 
the state \ipt) with respect to the control Ut- Here we choose the latter approach, where 
we have to solve the following equations 

a,|5V) +^^K)|5V) = - 2^5H{u°)\6^P) - z5'H{u°,M), 
d,\S4^)+iH{u°)m =-zSH{u°m), (6) 

which apply to any quantum closed system. Here \ijj°) is the solution of the Schrodinger 
equation for H{u°), and SH{u°), 5'^H{u°) are the Gateaux derivatives of the system 
Hamiltonian (defined as: 5H = — da° " U=o)- These derivatives are always 
analytically computable, since the dependence of the system Hamiltonian H{ut) on 
the control pulse Ut is always known, whereas the analytical dependence of both \6'^ip) 
and \6ijj) on Ut is known only in few cases (e.g., the driven and parametric harmonic 
oscillator |25]). Because of the latter, we need to solve the equations IQ which allow 
us to determine the matrix H. However, depending on the particular control problem, 
the BFGS method might be more efficient. Beside this, we note that the case of a 
linear quadratic regulator (or even Gaussian) control, that is, a system for which the 
state equation is linear and the performance criterion to be minimized is a quadratic 
form of the state and eventually also of the control pulse [19] , is not contemplated in 
our scheme. Indeed, the system state {ipt) is always dependent on the control pulse 
Ut through the Schrodinger equation of motion, even in the simple scenario where the 
system Hamiltonian H{ut) depends linearly on ut- Indeed, we are interested in the 
Hessian of the reduced cost functional J'[m] = ^[u, ip{u)], whose dependence on ut might 
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be not trivial through \ilj{ut)). Hence, the cases in which the Hessian rehes only upon 
the state do not concern J'[u]. 

We also remark that the above outlined formalism concerns state vectors in Hilbert 
spaces. If we would be interested in the optimization of unitary transformations 
U{t) = U{ut), then we should perform the Taylor expansion of F{U(T)) = itr{?7^?7(T)}, 
where f/g is the ideal unitary we wish to accomplish, and d is the dimension of the Hilbert 
space. Thus, the previuos analysis can be easily generalized to unitary operations. 

3. Applications of the method 

Let us first apply our method to an exactly solvable model. We consider the transport 
of a particle in a movable one-dimensional harmonic trap potential, for which the OC 
pulse and the functional dependence of the state on the controller are analytically known 
[22] . Such control problem is of relevance for the realization of a quantum ion processor. 
Indeed, optimistic estimates show that transport processes may account for 95% of the 
operation time of a quantum computation [2S]- Beside this, the harmonicity of the 
ion confinement in segmented Paul traps has been proven both numerically |2Z] and 
experimentally [28] . 



The system Hamiltonian is given by H{ut) = |{p^ + (x — UtY} with [x,p] 



I 



(we use harmonic oscillator units). We focus our attention on the ground state, that 
is, when the particle is initially prepared in the lowest vibrational state of the trap. 
The aim is to transport such a state over the distance Ax in a time T such that 
iPt{x) = e'^'^(j)o{x — Ax) = ipg{x), where (p is an unimportant phase factor and (f)o{x) is 
the Gaussian harmonic oscillator ground state wavefunction. Thus, our goal is to find a 
prescription such that when the OC pulse is perturbed the system has to reach at least 
the - a priori fixed - value of fidelity J-" G [0, 1]. 

In Ref. [22] the analytic solution for the time evolved state is provided. This 
enables us to compute analytically the Gateaux derivatives of the state of the system 
without the need of solving ([6]). In Fig. |2]^a) we show results for a distortion given 
by: 6ut = asm{K2Trt/T)u°, with u° being the OC pulse of Ref. [22] [see Eq. (5) 
therein]. Such simple distortion modulates the OC pulse at the rate k and gives a direct 
quantitative measure of the distortion degree applied to n°: the larger a is, the larger 
the infidelity. Thus, Fig. [2]^a) shows which is the largest admissible value of a for a 
fixed value of the (reduced) cost functional, whereas in the inset we show the deviation 
from the linear behaviour of 2A^J' = ^uH^u""" for k = 1. 

Now we write the cost functional as 2J'[ni] ~ ^uH^u"*". In this specific example it 
turns out, numerically, that the matrix elements of H are of the form Hnk = h + SHnk, 
with h,SHnk £ IR such that SHnk/h <^ 1. This means that the matrix elements are 
almost equal to each other. The eigenvalue problem for such a matrix can be well 
approximated by A^~^[A — {N — 2)h] = 0, that is, only an eigenvalue is non-zero. 
Thus, by defining h = Ylnk ^nk/{N — 2)^ we get for the non-vanishing eigenvalue 
the simple expression A0 = {N — 2)h. Given these remarks, the norm (jll) reduces 
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Figure 2. (Color online) (a) Distortion amplitude a (see text) vs. infidelity. Inset: 
second order Taylor expansion of J' {2A'^y[u°,Sut,Sut] = SuHSu^) vs. the reduced 
cost functional itself (red, k = 1) for different values of the distortion amplitude a, 
whereas the black line is a guide to the eye that helps to see when the second order 
approximation holds, (b) Infidelity tolerance parameter at', the solid (black) line 
represents [2X^{A'^ y [u°, 6ut, 6ut])'^ /Tr]^ , whereas the red (dashdot) line is a fit. 



to(J')~Ar0v/7r«3/(2A0). 

Finally, in order to determine at we proceed in the following way: since the norm 
is the "average cost functional" and at has to be unique, we simply use the inverted 
formula at = [2A0(J')^/(7rA/0)]3 for a given choice of Sut, and perform the substitution 
(J') — )■ (5uH(5u^/2 ~ 1 — F^i/jt)- (For several non-zero eigenvalues of H we would 
have had at = [2(J')^(X]fc=i V^-^fe/V^)^^]^-) Then, we choose some distortion 6ut 
and by varying the strength of such a distortion we collect the values of at versus the 
numerically exact overlap infidelities (i.e., without second order approximation), which 
is basically identified with J7 = 1 — J-", the fixed error threshold. We tested, however, 
that different kinds of distortions (e.g., Fourier-like model) produce practically the same 
curve as the one shown in Fig. Wljo)- This is easily understood, since what is relevant 
for the overlap infidelity terminal functional is the final state \iPt)- Hence, our method 
is general and it applies to any kind of distortion model. The reason for choosing the 
single frequency noise for the results displayed in Fig. M^a), was only to show the impact 
of the distortion on the cost functional in a simple and analytical way. 

Secondly, we perform a fit of the obtained curve for at as a function of the 
overlap infidelity. For the present example, showed in Fig. MJo), we found that the 
following function well represents the data: at = a(l — J-") + 6-\/l — J^, with a = 0.130 
and b = 0.029. Given that, all distortions of the optimal control pulse that satisfy 
the inequality ^uH^u""" < J nalM^ / {2\(ii) 



i{J^) for a fixed (a priori) value of 
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Figure 3. (Color online) Cost functional for the Landau-Zener model for 50 random 
realizations of the distortion Su{t) = asm{K2TTt/T)u°{t): the red (dashdot) line is the 
exact infidelity, the black (solid) line correspondes to the second order approximation 
I (see text), and the blue (thick) line is the infidelity threshold 1 — J^ — 0.01. 



desired fidelity J-", will yield an overlap fidelity f{ipT) > J^- From an operational 
point of view, this means that if we randomly choose a distortion 6ut such that 
I 



[duUdvJ ^/2\fD/{J\f^'ii)]'^^^/at < 1 — -7^, than we will certainly attain a state I^Jjt) 
whose overlap fidelity is greater than J-" (see also Fig. [T]). This permits us to find 
control pulses that might be more appropriate for an experimental implementation. 
This conclusion follows from the procedure we established for the determination of the 
infidelity tolerance Of We underscore, however, that the threshold i{J^) depends only 
on n°, because H relies only on u°, and that it is the result of an average in U through 
the path integral we defined. This provides a "universal" character to i{J^) for the given 
control problem. 

As a second example we consider the Landau-Zener model which has been proven 
useful to describe the tunneling of Bose-Einstein condensates in accelerated optical 
lattices [29] and the dynamics of a quench-induced phase transition in the quantum 
Ising model [30]. The model is described by the following system Hamiltonian: 
H{ut) = UtCFz + ^CTx, where cTz^cTx are Pauli matrices. The goal is to bring the system 
from the ground state \iI)q) of the Hamiltonian H{uq) to the ground state l^/^f.) of the 
Hamiltonian H{ut) through the avoided level crossing. As objective functional we 
consider y[ut\ = 1 — 'Re[{ip^\ipT)], that is, we also control the phase of the state. 
Even though the analytical dependence of the system Hamiltonian H{ut) is known, 
and therefore the Gateaux derivative {5H = dz^Ut), for such a control problem both the 
u° and the dependence of \ipt) on ut are not analytically known. Concerning the latter, 
this implies that is not possible to compute analytically the Gateaux derivatives \Sip) 
and \S'^ip), and therefore the Hessian H. Hence, we need to solve (^. Given that, we 
seek an u° by using the Krotov iterative method [H [31] . We then proceed on, as in the 
former example, by determining the infidelity tolerance at, for which we get a similar 
fit at = a(l -T)+ b{l - J^f with a = 0.016, b = 0.047, and c = 0.650. The criterion 
to be satisfied is then again given by: 5uH5u"'" < i{J^), but with a different numerical 
value for i{J^). 

In Fig. [3] it is showed the exact result of the infidelity (red) and the second order 
approximation X (black) for 50 random realizations of 6ut for a distortion similar to the 
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one of Fig. I21^a). Instead, the desired upper limit of infidelity 1 — J-" = 0.01, that the 
system has at least to attain, is represented by the horizontal blue line. Analogously 
to the former example, only the realizations of 5ut that fulfil (JuH^u"'" < £(0.99) have 
an infidelity below 1%. As for the previous example, the single frequency model has 
been adopted for convenience in order to illustrate the effectiveness of our approach as 
depicted in Fig. [31 Such a choice, however, does not invalidate our methodology, which 
we have tested for a Fourier-like distortion model (i.e. a sum over a finite number of 
harmonics with different amplitudes), similarly to the above outlined determination of 
the infidelity tolerance. 

In general, in order to identify robust control pulses amenable to an experimental 
implementation one should proceed as follows: 1) determine the optimal control 
u° with some numerical optimization algorithm; 2) compute the eigenvalues of the 
Hessian matrix H; 3) determine the infidelity tolerance at as described in the 
first example; 4) randomly generate distortions 5\i such that the condition X = 
[V2/(™?)(EfcliA/'fc/^A;)-^5uH5uT]2/3 < J is fulfilled for a fixed (a priori) value 
of J'] 5) define the new control pulse as u = u° + 6u; 6) choose the most suitable 
control pulse to be employed in the laboratory from the ensemble £ = {u}<j. We 
note, however, that such a "recipe" does not guarantee that we are able to obtain with 
certainty a control pulse more amenable for use in an experiment, but at least it helps 
to design new ones that are still able to yield a value of the cost functional below the 
fixed threshold J'. In other words, the set <£^ fl ^ might be empty, and the elements of 
£ are not necessarily close to the optimal u° (e.g., see Ref. [22])- In this respect, such 
a procedure might help the quest of both robust and experimentally feasible pulses for 
controlling different quantum phenomena. 

4. Conclusions 

In conclusion we have presented a method to assess the error of solutions to OC 
problems. The method might be a helpful tool for experimentalists in order to design 
experiments robust against source of noise. For instance, to estimate how much the 
schemes for realizing quantum gates are robust against imperfections of the optimal pulse 
shape. Compared to other methods, such as the one of Ref. [13], our technique does not 
need the simulation of a large ensemble of samples in order to minimize on average both 
the mean and the variance of the objective functional. Instead, once the Hessian matrix 
H and the infidelity tolerance at are known, which only rely on the optimal (known) 
control pulses, one has simply to randomly generate the distortion 6u and perform the 
matrix multiplication 5uH5u^, which is a less demanding computational task than the 
application of a genetic algorithm, as the one proposed in Ref. pJJ. Besides this, our 
method can be applied to both known error models and to evaluate unknown errors 
such as random telegraphic noise. For the future, we plan to extend our formalism 
to dissipative quantum systems (e.g., systems governed by the Born-Markov master 
equation) . 
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